home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dhgeqz.f < prev    next >
Text File  |  1996-11-04  |  42KB  |  1,234 lines
  1.       SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB,
  2.      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
  3.      $                   LWORK, INFO )
  4. *
  5. *  -- LAPACK routine (version 2.0) --
  6. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  7. *     Courant Institute, Argonne National Lab, and Rice University
  8. *     September 30, 1994
  9. *
  10. *     .. Scalar Arguments ..
  11.       CHARACTER          COMPQ, COMPZ, JOB
  12.       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
  13. *     ..
  14. *     .. Array Arguments ..
  15.       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
  16.      $                   B( LDB, * ), BETA( * ), Q( LDQ, * ), WORK( * ),
  17.      $                   Z( LDZ, * )
  18. *     ..
  19. *
  20. *  Purpose
  21. *  =======
  22. *
  23. *  DHGEQZ implements a single-/double-shift version of the QZ method for
  24. *  finding the generalized eigenvalues
  25. *
  26. *  w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j)   of the equation
  27. *
  28. *       det( A - w(i) B ) = 0
  29. *
  30. *  In addition, the pair A,B may be reduced to generalized Schur form:
  31. *  B is upper triangular, and A is block upper triangular, where the
  32. *  diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having
  33. *  complex generalized eigenvalues (see the description of the argument
  34. *  JOB.)
  35. *
  36. *  If JOB='S', then the pair (A,B) is simultaneously reduced to Schur
  37. *  form by applying one orthogonal tranformation (usually called Q) on
  38. *  the left and another (usually called Z) on the right.  The 2-by-2
  39. *  upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks
  40. *  of A will be reduced to positive diagonal matrices.  (I.e.,
  41. *  if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and
  42. *  B(j+1,j+1) will be positive.)
  43. *
  44. *  If JOB='E', then at each iteration, the same transformations
  45. *  are computed, but they are only applied to those parts of A and B
  46. *  which are needed to compute ALPHAR, ALPHAI, and BETAR.
  47. *
  48. *  If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal
  49. *  transformations used to reduce (A,B) are accumulated into the arrays
  50. *  Q and Z s.t.:
  51. *
  52. *       Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
  53. *       Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*
  54. *
  55. *  Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
  56. *       Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
  57. *       pp. 241--256.
  58. *
  59. *  Arguments
  60. *  =========
  61. *
  62. *  JOB     (input) CHARACTER*1
  63. *          = 'E': compute only ALPHAR, ALPHAI, and BETA.  A and B will
  64. *                 not necessarily be put into generalized Schur form.
  65. *          = 'S': put A and B into generalized Schur form, as well
  66. *                 as computing ALPHAR, ALPHAI, and BETA.
  67. *
  68. *  COMPQ   (input) CHARACTER*1
  69. *          = 'N': do not modify Q.
  70. *          = 'V': multiply the array Q on the right by the transpose of
  71. *                 the orthogonal tranformation that is applied to the
  72. *                 left side of A and B to reduce them to Schur form.
  73. *          = 'I': like COMPQ='V', except that Q will be initialized to
  74. *                 the identity first.
  75. *
  76. *  COMPZ   (input) CHARACTER*1
  77. *          = 'N': do not modify Z.
  78. *          = 'V': multiply the array Z on the right by the orthogonal
  79. *                 tranformation that is applied to the right side of
  80. *                 A and B to reduce them to Schur form.
  81. *          = 'I': like COMPZ='V', except that Z will be initialized to
  82. *                 the identity first.
  83. *
  84. *  N       (input) INTEGER
  85. *          The order of the matrices A, B, Q, and Z.  N >= 0.
  86. *
  87. *  ILO     (input) INTEGER
  88. *  IHI     (input) INTEGER
  89. *          It is assumed that A is already upper triangular in rows and
  90. *          columns 1:ILO-1 and IHI+1:N.
  91. *          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
  92. *
  93. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
  94. *          On entry, the N-by-N upper Hessenberg matrix A.  Elements
  95. *          below the subdiagonal must be zero.
  96. *          If JOB='S', then on exit A and B will have been
  97. *             simultaneously reduced to generalized Schur form.
  98. *          If JOB='E', then on exit A will have been destroyed.
  99. *             The diagonal blocks will be correct, but the off-diagonal
  100. *             portion will be meaningless.
  101. *
  102. *  LDA     (input) INTEGER
  103. *          The leading dimension of the array A.  LDA >= max( 1, N ).
  104. *
  105. *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
  106. *          On entry, the N-by-N upper triangular matrix B.  Elements
  107. *          below the diagonal must be zero.  2-by-2 blocks in B
  108. *          corresponding to 2-by-2 blocks in A will be reduced to
  109. *          positive diagonal form.  (I.e., if A(j+1,j) is non-zero,
  110. *          then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be
  111. *          positive.)
  112. *          If JOB='S', then on exit A and B will have been
  113. *             simultaneously reduced to Schur form.
  114. *          If JOB='E', then on exit B will have been destroyed.
  115. *             Elements corresponding to diagonal blocks of A will be
  116. *             correct, but the off-diagonal portion will be meaningless.
  117. *
  118. *  LDB     (input) INTEGER
  119. *          The leading dimension of the array B.  LDB >= max( 1, N ).
  120. *
  121. *  ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
  122. *          ALPHAR(1:N) will be set to real parts of the diagonal
  123. *          elements of A that would result from reducing A and B to
  124. *          Schur form and then further reducing them both to triangular
  125. *          form using unitary transformations s.t. the diagonal of B
  126. *          was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block
  127. *          (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j).
  128. *          Note that the (real or complex) values
  129. *          (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
  130. *          generalized eigenvalues of the matrix pencil A - wB.
  131. *
  132. *  ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
  133. *          ALPHAI(1:N) will be set to imaginary parts of the diagonal
  134. *          elements of A that would result from reducing A and B to
  135. *          Schur form and then further reducing them both to triangular
  136. *          form using unitary transformations s.t. the diagonal of B
  137. *          was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block
  138. *          (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0.
  139. *          Note that the (real or complex) values
  140. *          (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
  141. *          generalized eigenvalues of the matrix pencil A - wB.
  142. *
  143. *  BETA    (output) DOUBLE PRECISION array, dimension (N)
  144. *          BETA(1:N) will be set to the (real) diagonal elements of B
  145. *          that would result from reducing A and B to Schur form and
  146. *          then further reducing them both to triangular form using
  147. *          unitary transformations s.t. the diagonal of B was
  148. *          non-negative real.  Thus, if A(j,j) is in a 1-by-1 block
  149. *          (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j).
  150. *          Note that the (real or complex) values
  151. *          (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
  152. *          generalized eigenvalues of the matrix pencil A - wB.
  153. *          (Note that BETA(1:N) will always be non-negative, and no
  154. *          BETAI is necessary.)
  155. *
  156. *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
  157. *          If COMPQ='N', then Q will not be referenced.
  158. *          If COMPQ='V' or 'I', then the transpose of the orthogonal
  159. *             transformations which are applied to A and B on the left
  160. *             will be applied to the array Q on the right.
  161. *
  162. *  LDQ     (input) INTEGER
  163. *          The leading dimension of the array Q.  LDQ >= 1.
  164. *          If COMPQ='V' or 'I', then LDQ >= N.
  165. *
  166. *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
  167. *          If COMPZ='N', then Z will not be referenced.
  168. *          If COMPZ='V' or 'I', then the orthogonal transformations
  169. *             which are applied to A and B on the right will be applied
  170. *             to the array Z on the right.
  171. *
  172. *  LDZ     (input) INTEGER
  173. *          The leading dimension of the array Z.  LDZ >= 1.
  174. *          If COMPZ='V' or 'I', then LDZ >= N.
  175. *
  176. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  177. *          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
  178. *
  179. *  LWORK   (input) INTEGER
  180. *          The dimension of the array WORK.  LWORK >= max(1,N).
  181. *
  182. *  INFO    (output) INTEGER
  183. *          = 0: successful exit
  184. *          < 0: if INFO = -i, the i-th argument had an illegal value
  185. *          = 1,...,N: the QZ iteration did not converge.  (A,B) is not
  186. *                     in Schur form, but ALPHAR(i), ALPHAI(i), and
  187. *                     BETA(i), i=INFO+1,...,N should be correct.
  188. *          = N+1,...,2*N: the shift calculation failed.  (A,B) is not
  189. *                     in Schur form, but ALPHAR(i), ALPHAI(i), and
  190. *                     BETA(i), i=INFO-N+1,...,N should be correct.
  191. *          > 2*N:     various "impossible" errors.
  192. *
  193. *  Further Details
  194. *  ===============
  195. *
  196. *  Iteration counters:
  197. *
  198. *  JITER  -- counts iterations.
  199. *  IITER  -- counts iterations run since ILAST was last
  200. *            changed.  This is therefore reset only when a 1-by-1 or
  201. *            2-by-2 block deflates off the bottom.
  202. *
  203. *  =====================================================================
  204. *
  205. *     .. Parameters ..
  206. *    $                     SAFETY = 1.0E+0 )
  207.       DOUBLE PRECISION   HALF, ZERO, ONE, SAFETY
  208.       PARAMETER          ( HALF = 0.5D+0, ZERO = 0.0D+0, ONE = 1.0D+0,
  209.      $                   SAFETY = 1.0D+2 )
  210. *     ..
  211. *     .. Local Scalars ..
  212.       LOGICAL            ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ
  213.       INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
  214.      $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
  215.      $                   JR, MAXIT
  216.       DOUBLE PRECISION   A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
  217.      $                   AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
  218.      $                   AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
  219.      $                   B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
  220.      $                   BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
  221.      $                   CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
  222.      $                   SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T,
  223.      $                   TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
  224.      $                   U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
  225.      $                   WR2
  226. *     ..
  227. *     .. Local Arrays ..
  228.       DOUBLE PRECISION   V( 3 )
  229. *     ..
  230. *     .. External Functions ..
  231.       LOGICAL            LSAME
  232.       DOUBLE PRECISION   DLAMCH, DLANHS, DLAPY2, DLAPY3
  233.       EXTERNAL           LSAME, DLAMCH, DLANHS, DLAPY2, DLAPY3
  234. *     ..
  235. *     .. External Subroutines ..
  236.       EXTERNAL           DLAG2, DLARFG, DLARTG, DLASET, DLASV2, DROT,
  237.      $                   XERBLA
  238. *     ..
  239. *     .. Intrinsic Functions ..
  240.       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
  241. *     ..
  242. *     .. Executable Statements ..
  243. *
  244. *     Decode JOB, COMPQ, COMPZ
  245. *
  246.       IF( LSAME( JOB, 'E' ) ) THEN
  247.          ILSCHR = .FALSE.
  248.          ISCHUR = 1
  249.       ELSE IF( LSAME( JOB, 'S' ) ) THEN
  250.          ILSCHR = .TRUE.
  251.          ISCHUR = 2
  252.       ELSE
  253.          ISCHUR = 0
  254.       END IF
  255. *
  256.       IF( LSAME( COMPQ, 'N' ) ) THEN
  257.          ILQ = .FALSE.
  258.          ICOMPQ = 1
  259.       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
  260.          ILQ = .TRUE.
  261.          ICOMPQ = 2
  262.       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
  263.          ILQ = .TRUE.
  264.          ICOMPQ = 3
  265.       ELSE
  266.          ICOMPQ = 0
  267.       END IF
  268. *
  269.       IF( LSAME( COMPZ, 'N' ) ) THEN
  270.          ILZ = .FALSE.
  271.          ICOMPZ = 1
  272.       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
  273.          ILZ = .TRUE.
  274.          ICOMPZ = 2
  275.       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
  276.          ILZ = .TRUE.
  277.          ICOMPZ = 3
  278.       ELSE
  279.          ICOMPZ = 0
  280.       END IF
  281. *
  282. *     Check Argument Values
  283. *
  284.       INFO = 0
  285.       IF( ISCHUR.EQ.0 ) THEN
  286.          INFO = -1
  287.       ELSE IF( ICOMPQ.EQ.0 ) THEN
  288.          INFO = -2
  289.       ELSE IF( ICOMPZ.EQ.0 ) THEN
  290.          INFO = -3
  291.       ELSE IF( N.LT.0 ) THEN
  292.          INFO = -4
  293.       ELSE IF( ILO.LT.1 ) THEN
  294.          INFO = -5
  295.       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
  296.          INFO = -6
  297.       ELSE IF( LDA.LT.N ) THEN
  298.          INFO = -8
  299.       ELSE IF( LDB.LT.N ) THEN
  300.          INFO = -10
  301.       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
  302.          INFO = -15
  303.       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
  304.          INFO = -17
  305.       ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN
  306.          INFO = -19
  307.       END IF
  308.       IF( INFO.NE.0 ) THEN
  309.          CALL XERBLA( 'DHGEQZ', -INFO )
  310.          RETURN
  311.       END IF
  312. *
  313. *     Quick return if possible
  314. *
  315.       IF( N.LE.0 ) THEN
  316.          WORK( 1 ) = DBLE( 1 )
  317.          RETURN
  318.       END IF
  319. *
  320. *     Initialize Q and Z
  321. *
  322.       IF( ICOMPQ.EQ.3 )
  323.      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
  324.       IF( ICOMPZ.EQ.3 )
  325.      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
  326. *
  327. *     Machine Constants
  328. *
  329.       IN = IHI + 1 - ILO
  330.       SAFMIN = DLAMCH( 'S' )
  331.       SAFMAX = ONE / SAFMIN
  332.       ULP = DLAMCH( 'E' )*DLAMCH( 'B' )
  333.       ANORM = DLANHS( 'F', IN, A( ILO, ILO ), LDA, WORK )
  334.       BNORM = DLANHS( 'F', IN, B( ILO, ILO ), LDB, WORK )
  335.       ATOL = MAX( SAFMIN, ULP*ANORM )
  336.       BTOL = MAX( SAFMIN, ULP*BNORM )
  337.       ASCALE = ONE / MAX( SAFMIN, ANORM )
  338.       BSCALE = ONE / MAX( SAFMIN, BNORM )
  339. *
  340. *     Set Eigenvalues IHI+1:N
  341. *
  342.       DO 30 J = IHI + 1, N
  343.          IF( B( J, J ).LT.ZERO ) THEN
  344.             IF( ILSCHR ) THEN
  345.                DO 10 JR = 1, J
  346.                   A( JR, J ) = -A( JR, J )
  347.                   B( JR, J ) = -B( JR, J )
  348.    10          CONTINUE
  349.             ELSE
  350.                A( J, J ) = -A( J, J )
  351.                B( J, J ) = -B( J, J )
  352.             END IF
  353.             IF( ILZ ) THEN
  354.                DO 20 JR = 1, N
  355.                   Z( JR, J ) = -Z( JR, J )
  356.    20          CONTINUE
  357.             END IF
  358.          END IF
  359.          ALPHAR( J ) = A( J, J )
  360.          ALPHAI( J ) = ZERO
  361.          BETA( J ) = B( J, J )
  362.    30 CONTINUE
  363. *
  364. *     If IHI < ILO, skip QZ steps
  365. *
  366.       IF( IHI.LT.ILO )
  367.      $   GO TO 380
  368. *
  369. *     MAIN QZ ITERATION LOOP
  370. *
  371. *     Initialize dynamic indices
  372. *
  373. *     Eigenvalues ILAST+1:N have been found.
  374. *        Column operations modify rows IFRSTM:whatever.
  375. *        Row operations modify columns whatever:ILASTM.
  376. *
  377. *     If only eigenvalues are being computed, then
  378. *        IFRSTM is the row of the last splitting row above row ILAST;
  379. *        this is always at least ILO.
  380. *     IITER counts iterations since the last eigenvalue was found,
  381. *        to tell when to use an extraordinary shift.
  382. *     MAXIT is the maximum number of QZ sweeps allowed.
  383. *
  384.       ILAST = IHI
  385.       IF( ILSCHR ) THEN
  386.          IFRSTM = 1
  387.          ILASTM = N
  388.       ELSE
  389.          IFRSTM = ILO
  390.          ILASTM = IHI
  391.       END IF
  392.       IITER = 0
  393.       ESHIFT = ZERO
  394.       MAXIT = 30*( IHI-ILO+1 )
  395. *
  396.       DO 360 JITER = 1, MAXIT
  397. *
  398. *        Split the matrix if possible.
  399. *
  400. *        Two tests:
  401. *           1: A(j,j-1)=0  or  j=ILO
  402. *           2: B(j,j)=0
  403. *
  404.          IF( ILAST.EQ.ILO ) THEN
  405. *
  406. *           Special case: j=ILAST
  407. *
  408.             GO TO 80
  409.          ELSE
  410.             IF( ABS( A( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
  411.                A( ILAST, ILAST-1 ) = ZERO
  412.                GO TO 80
  413.             END IF
  414.          END IF
  415. *
  416.          IF( ABS( B( ILAST, ILAST ) ).LE.BTOL ) THEN
  417.             B( ILAST, ILAST ) = ZERO
  418.             GO TO 70
  419.          END IF
  420. *
  421. *        General case: j<ILAST
  422. *
  423.          DO 60 J = ILAST - 1, ILO, -1
  424. *
  425. *           Test 1: for A(j,j-1)=0 or j=ILO
  426. *
  427.             IF( J.EQ.ILO ) THEN
  428.                ILAZRO = .TRUE.
  429.             ELSE
  430.                IF( ABS( A( J, J-1 ) ).LE.ATOL ) THEN
  431.                   A( J, J-1 ) = ZERO
  432.                   ILAZRO = .TRUE.
  433.                ELSE
  434.                   ILAZRO = .FALSE.
  435.                END IF
  436.             END IF
  437. *
  438. *           Test 2: for B(j,j)=0
  439. *
  440.             IF( ABS( B( J, J ) ).LT.BTOL ) THEN
  441.                B( J, J ) = ZERO
  442. *
  443. *              Test 1a: Check for 2 consecutive small subdiagonals in A
  444. *
  445.                ILAZR2 = .FALSE.
  446.                IF( .NOT.ILAZRO ) THEN
  447.                   TEMP = ABS( A( J, J-1 ) )
  448.                   TEMP2 = ABS( A( J, J ) )
  449.                   TEMPR = MAX( TEMP, TEMP2 )
  450.                   IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
  451.                      TEMP = TEMP / TEMPR
  452.                      TEMP2 = TEMP2 / TEMPR
  453.                   END IF
  454.                   IF( TEMP*( ASCALE*ABS( A( J+1, J ) ) ).LE.TEMP2*
  455.      $                ( ASCALE*ATOL ) )ILAZR2 = .TRUE.
  456.                END IF
  457. *
  458. *              If both tests pass (1 & 2), i.e., the leading diagonal
  459. *              element of B in the block is zero, split a 1x1 block off
  460. *              at the top. (I.e., at the J-th row/column) The leading
  461. *              diagonal element of the remainder can also be zero, so
  462. *              this may have to be done repeatedly.
  463. *
  464.                IF( ILAZRO .OR. ILAZR2 ) THEN
  465.                   DO 40 JCH = J, ILAST - 1
  466.                      TEMP = A( JCH, JCH )
  467.                      CALL DLARTG( TEMP, A( JCH+1, JCH ), C, S,
  468.      $                            A( JCH, JCH ) )
  469.                      A( JCH+1, JCH ) = ZERO
  470.                      CALL DROT( ILASTM-JCH, A( JCH, JCH+1 ), LDA,
  471.      $                          A( JCH+1, JCH+1 ), LDA, C, S )
  472.                      CALL DROT( ILASTM-JCH, B( JCH, JCH+1 ), LDB,
  473.      $                          B( JCH+1, JCH+1 ), LDB, C, S )
  474.                      IF( ILQ )
  475.      $                  CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
  476.      $                             C, S )
  477.                      IF( ILAZR2 )
  478.      $                  A( JCH, JCH-1 ) = A( JCH, JCH-1 )*C
  479.                      ILAZR2 = .FALSE.
  480.                      IF( ABS( B( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
  481.                         IF( JCH+1.GE.ILAST ) THEN
  482.                            GO TO 80
  483.                         ELSE
  484.                            IFIRST = JCH + 1
  485.                            GO TO 110
  486.                         END IF
  487.                      END IF
  488.                      B( JCH+1, JCH+1 ) = ZERO
  489.    40             CONTINUE
  490.                   GO TO 70
  491.                ELSE
  492. *
  493. *                 Only test 2 passed -- chase the zero to B(ILAST,ILAST)
  494. *                 Then process as in the case B(ILAST,ILAST)=0
  495. *
  496.                   DO 50 JCH = J, ILAST - 1
  497.                      TEMP = B( JCH, JCH+1 )
  498.                      CALL DLARTG( TEMP, B( JCH+1, JCH+1 ), C, S,
  499.      $                            B( JCH, JCH+1 ) )
  500.                      B( JCH+1, JCH+1 ) = ZERO
  501.                      IF( JCH.LT.ILASTM-1 )
  502.      $                  CALL DROT( ILASTM-JCH-1, B( JCH, JCH+2 ), LDB,
  503.      $                             B( JCH+1, JCH+2 ), LDB, C, S )
  504.                      CALL DROT( ILASTM-JCH+2, A( JCH, JCH-1 ), LDA,
  505.      $                          A( JCH+1, JCH-1 ), LDA, C, S )
  506.                      IF( ILQ )
  507.      $                  CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
  508.      $                             C, S )
  509.                      TEMP = A( JCH+1, JCH )
  510.                      CALL DLARTG( TEMP, A( JCH+1, JCH-1 ), C, S,
  511.      $                            A( JCH+1, JCH ) )
  512.                      A( JCH+1, JCH-1 ) = ZERO
  513.                      CALL DROT( JCH+1-IFRSTM, A( IFRSTM, JCH ), 1,
  514.      $                          A( IFRSTM, JCH-1 ), 1, C, S )
  515.                      CALL DROT( JCH-IFRSTM, B( IFRSTM, JCH ), 1,
  516.      $                          B( IFRSTM, JCH-1 ), 1, C, S )
  517.                      IF( ILZ )
  518.      $                  CALL DROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
  519.      $                             C, S )
  520.    50             CONTINUE
  521.                   GO TO 70
  522.                END IF
  523.             ELSE IF( ILAZRO ) THEN
  524. *
  525. *              Only test 1 passed -- work on J:ILAST
  526. *
  527.                IFIRST = J
  528.                GO TO 110
  529.             END IF
  530. *
  531. *           Neither test passed -- try next J
  532. *
  533.    60    CONTINUE
  534. *
  535. *        (Drop-through is "impossible")
  536. *
  537.          INFO = N + 1
  538.          GO TO 420
  539. *
  540. *        B(ILAST,ILAST)=0 -- clear A(ILAST,ILAST-1) to split off a
  541. *        1x1 block.
  542. *
  543.    70    CONTINUE
  544.          TEMP = A( ILAST, ILAST )
  545.          CALL DLARTG( TEMP, A( ILAST, ILAST-1 ), C, S,
  546.      $                A( ILAST, ILAST ) )
  547.          A( ILAST, ILAST-1 ) = ZERO
  548.          CALL DROT( ILAST-IFRSTM, A( IFRSTM, ILAST ), 1,
  549.      $              A( IFRSTM, ILAST-1 ), 1, C, S )
  550.          CALL DROT( ILAST-IFRSTM, B( IFRSTM, ILAST ), 1,
  551.      $              B( IFRSTM, ILAST-1 ), 1, C, S )
  552.          IF( ILZ )
  553.      $      CALL DROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
  554. *
  555. *        A(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
  556. *                              and BETA
  557. *
  558.    80    CONTINUE
  559.          IF( B( ILAST, ILAST ).LT.ZERO ) THEN
  560.             IF( ILSCHR ) THEN
  561.                DO 90 J = IFRSTM, ILAST
  562.                   A( J, ILAST ) = -A( J, ILAST )
  563.                   B( J, ILAST ) = -B( J, ILAST )
  564.    90          CONTINUE
  565.             ELSE
  566.                A( ILAST, ILAST ) = -A( ILAST, ILAST )
  567.                B( ILAST, ILAST ) = -B( ILAST, ILAST )
  568.             END IF
  569.             IF( ILZ ) THEN
  570.                DO 100 J = 1, N
  571.                   Z( J, ILAST ) = -Z( J, ILAST )
  572.   100          CONTINUE
  573.             END IF
  574.          END IF
  575.          ALPHAR( ILAST ) = A( ILAST, ILAST )
  576.          ALPHAI( ILAST ) = ZERO
  577.          BETA( ILAST ) = B( ILAST, ILAST )
  578. *
  579. *        Go to next block -- exit if finished.
  580. *
  581.          ILAST = ILAST - 1
  582.          IF( ILAST.LT.ILO )
  583.      $      GO TO 380
  584. *
  585. *        Reset counters
  586. *
  587.          IITER = 0
  588.          ESHIFT = ZERO
  589.          IF( .NOT.ILSCHR ) THEN
  590.             ILASTM = ILAST
  591.             IF( IFRSTM.GT.ILAST )
  592.      $         IFRSTM = ILO
  593.          END IF
  594.          GO TO 350
  595. *
  596. *        QZ step
  597. *
  598. *        This iteration only involves rows/columns IFIRST:ILAST. We
  599. *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
  600. *
  601.   110    CONTINUE
  602.          IITER = IITER + 1
  603.          IF( .NOT.ILSCHR ) THEN
  604.             IFRSTM = IFIRST
  605.          END IF
  606. *
  607. *        Compute single shifts.
  608. *
  609. *        At this point, IFIRST < ILAST, and the diagonal elements of
  610. *        B(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
  611. *        magnitude)
  612. *
  613.          IF( ( IITER / 10 )*10.EQ.IITER ) THEN
  614. *
  615. *           Exceptional shift.  Chosen for no particularly good reason.
  616. *           (Single shift only.)
  617. *
  618.             IF( ( DBLE( MAXIT )*SAFMIN )*ABS( A( ILAST-1, ILAST ) ).LT.
  619.      $          ABS( B( ILAST-1, ILAST-1 ) ) ) THEN
  620.                ESHIFT = ESHIFT + A( ILAST-1, ILAST ) /
  621.      $                  B( ILAST-1, ILAST-1 )
  622.             ELSE
  623.                ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
  624.             END IF
  625.             S1 = ONE
  626.             WR = ESHIFT
  627. *
  628.          ELSE
  629. *
  630. *           Shifts based on the generalized eigenvalues of the
  631. *           bottom-right 2x2 block of A and B. The first eigenvalue
  632. *           returned by DLAG2 is the Wilkinson shift (AEP p.512),
  633. *
  634.             CALL DLAG2( A( ILAST-1, ILAST-1 ), LDA,
  635.      $                  B( ILAST-1, ILAST-1 ), LDB, SAFMIN*SAFETY, S1,
  636.      $                  S2, WR, WR2, WI )
  637. *
  638.             TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
  639.             IF( WI.NE.ZERO )
  640.      $         GO TO 200
  641.          END IF
  642. *
  643. *        Fiddle with shift to avoid overflow
  644. *
  645.          TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX )
  646.          IF( S1.GT.TEMP ) THEN
  647.             SCALE = TEMP / S1
  648.          ELSE
  649.             SCALE = ONE
  650.          END IF
  651. *
  652.          TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX )
  653.          IF( ABS( WR ).GT.TEMP )
  654.      $      SCALE = MIN( SCALE, TEMP / ABS( WR ) )
  655.          S1 = SCALE*S1
  656.          WR = SCALE*WR
  657. *
  658. *        Now check for two consecutive small subdiagonals.
  659. *
  660.          DO 120 J = ILAST - 1, IFIRST + 1, -1
  661.             ISTART = J
  662.             TEMP = ABS( S1*A( J, J-1 ) )
  663.             TEMP2 = ABS( S1*A( J, J )-WR*B( J, J ) )
  664.             TEMPR = MAX( TEMP, TEMP2 )
  665.             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
  666.                TEMP = TEMP / TEMPR
  667.                TEMP2 = TEMP2 / TEMPR
  668.             END IF
  669.             IF( ABS( ( ASCALE*A( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )*
  670.      $          TEMP2 )GO TO 130
  671.   120    CONTINUE
  672. *
  673.          ISTART = IFIRST
  674.   130    CONTINUE
  675. *
  676. *        Do an implicit single-shift QZ sweep.
  677. *
  678. *        Initial Q
  679. *
  680.          TEMP = S1*A( ISTART, ISTART ) - WR*B( ISTART, ISTART )
  681.          TEMP2 = S1*A( ISTART+1, ISTART )
  682.          CALL DLARTG( TEMP, TEMP2, C, S, TEMPR )
  683. *
  684. *        Sweep
  685. *
  686.          DO 190 J = ISTART, ILAST - 1
  687.             IF( J.GT.ISTART ) THEN
  688.                TEMP = A( J, J-1 )
  689.                CALL DLARTG( TEMP, A( J+1, J-1 ), C, S, A( J, J-1 ) )
  690.                A( J+1, J-1 ) = ZERO
  691.             END IF
  692. *
  693.             DO 140 JC = J, ILASTM
  694.                TEMP = C*A( J, JC ) + S*A( J+1, JC )
  695.                A( J+1, JC ) = -S*A( J, JC ) + C*A( J+1, JC )
  696.                A( J, JC ) = TEMP
  697.                TEMP2 = C*B( J, JC ) + S*B( J+1, JC )
  698.                B( J+1, JC ) = -S*B( J, JC ) + C*B( J+1, JC )
  699.                B( J, JC ) = TEMP2
  700.   140       CONTINUE
  701.             IF( ILQ ) THEN
  702.                DO 150 JR = 1, N
  703.                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
  704.                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
  705.                   Q( JR, J ) = TEMP
  706.   150          CONTINUE
  707.             END IF
  708. *
  709.             TEMP = B( J+1, J+1 )
  710.             CALL DLARTG( TEMP, B( J+1, J ), C, S, B( J+1, J+1 ) )
  711.             B( J+1, J ) = ZERO
  712. *
  713.             DO 160 JR = IFRSTM, MIN( J+2, ILAST )
  714.                TEMP = C*A( JR, J+1 ) + S*A( JR, J )
  715.                A( JR, J ) = -S*A( JR, J+1 ) + C*A( JR, J )
  716.                A( JR, J+1 ) = TEMP
  717.   160       CONTINUE
  718.             DO 170 JR = IFRSTM, J
  719.                TEMP = C*B( JR, J+1 ) + S*B( JR, J )
  720.                B( JR, J ) = -S*B( JR, J+1 ) + C*B( JR, J )
  721.                B( JR, J+1 ) = TEMP
  722.   170       CONTINUE
  723.             IF( ILZ ) THEN
  724.                DO 180 JR = 1, N
  725.                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
  726.                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
  727.                   Z( JR, J+1 ) = TEMP
  728.   180          CONTINUE
  729.             END IF
  730.   190    CONTINUE
  731. *
  732.          GO TO 350
  733. *
  734. *        Use Francis double-shift
  735. *
  736. *        Note: the Francis double-shift should work with real shifts,
  737. *              but only if the block is at least 3x3.
  738. *              This code may break if this point is reached with
  739. *              a 2x2 block with real eigenvalues.
  740. *
  741.   200    CONTINUE
  742.          IF( IFIRST+1.EQ.ILAST ) THEN
  743. *
  744. *           Special case -- 2x2 block with complex eigenvectors
  745. *
  746. *           Step 1: Standardize, that is, rotate so that
  747. *
  748. *                       ( B11  0  )
  749. *                   B = (         )  with B11 non-negative.
  750. *                       (  0  B22 )
  751. *
  752.             CALL DLASV2( B( ILAST-1, ILAST-1 ), B( ILAST-1, ILAST ),
  753.      $                   B( ILAST, ILAST ), B22, B11, SR, CR, SL, CL )
  754. *
  755.             IF( B11.LT.ZERO ) THEN
  756.                CR = -CR
  757.                SR = -SR
  758.                B11 = -B11
  759.                B22 = -B22
  760.             END IF
  761. *
  762.             CALL DROT( ILASTM+1-IFIRST, A( ILAST-1, ILAST-1 ), LDA,
  763.      $                 A( ILAST, ILAST-1 ), LDA, CL, SL )
  764.             CALL DROT( ILAST+1-IFRSTM, A( IFRSTM, ILAST-1 ), 1,
  765.      $                 A( IFRSTM, ILAST ), 1, CR, SR )
  766. *
  767.             IF( ILAST.LT.ILASTM )
  768.      $         CALL DROT( ILASTM-ILAST, B( ILAST-1, ILAST+1 ), LDB,
  769.      $                    B( ILAST, ILAST+1 ), LDA, CL, SL )
  770.             IF( IFRSTM.LT.ILAST-1 )
  771.      $         CALL DROT( IFIRST-IFRSTM, B( IFRSTM, ILAST-1 ), 1,
  772.      $                    B( IFRSTM, ILAST ), 1, CR, SR )
  773. *
  774.             IF( ILQ )
  775.      $         CALL DROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
  776.      $                    SL )
  777.             IF( ILZ )
  778.      $         CALL DROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR,
  779.      $                    SR )
  780. *
  781.             B( ILAST-1, ILAST-1 ) = B11
  782.             B( ILAST-1, ILAST ) = ZERO
  783.             B( ILAST, ILAST-1 ) = ZERO
  784.             B( ILAST, ILAST ) = B22
  785. *
  786. *           If B22 is negative, negate column ILAST
  787. *
  788.             IF( B22.LT.ZERO ) THEN
  789.                DO 210 J = IFRSTM, ILAST
  790.                   A( J, ILAST ) = -A( J, ILAST )
  791.                   B( J, ILAST ) = -B( J, ILAST )
  792.   210          CONTINUE
  793. *
  794.                IF( ILZ ) THEN
  795.                   DO 220 J = 1, N
  796.                      Z( J, ILAST ) = -Z( J, ILAST )
  797.   220             CONTINUE
  798.                END IF
  799.             END IF
  800. *
  801. *           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
  802. *
  803. *           Recompute shift
  804. *
  805.             CALL DLAG2( A( ILAST-1, ILAST-1 ), LDA,
  806.      $                  B( ILAST-1, ILAST-1 ), LDB, SAFMIN*SAFETY, S1,
  807.      $                  TEMP, WR, TEMP2, WI )
  808. *
  809. *           If standardization has perturbed the shift onto real line,
  810. *           do another (real single-shift) QR step.
  811. *
  812.             IF( WI.EQ.ZERO )
  813.      $         GO TO 350
  814.             S1INV = ONE / S1
  815. *
  816. *           Do EISPACK (QZVAL) computation of alpha and beta
  817. *
  818.             A11 = A( ILAST-1, ILAST-1 )
  819.             A21 = A( ILAST, ILAST-1 )
  820.             A12 = A( ILAST-1, ILAST )
  821.             A22 = A( ILAST, ILAST )
  822. *
  823. *           Compute complex Givens rotation on right
  824. *           (Assume some element of C = (sA - wB) > unfl )
  825. *                            __
  826. *           (sA - wB) ( CZ   -SZ )
  827. *                     ( SZ    CZ )
  828. *
  829.             C11R = S1*A11 - WR*B11
  830.             C11I = -WI*B11
  831.             C12 = S1*A12
  832.             C21 = S1*A21
  833.             C22R = S1*A22 - WR*B22
  834.             C22I = -WI*B22
  835. *
  836.             IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+
  837.      $          ABS( C22R )+ABS( C22I ) ) THEN
  838.                T = DLAPY3( C12, C11R, C11I )
  839.                CZ = C12 / T
  840.                SZR = -C11R / T
  841.                SZI = -C11I / T
  842.             ELSE
  843.                CZ = DLAPY2( C22R, C22I )
  844.                IF( CZ.LE.SAFMIN ) THEN
  845.                   CZ = ZERO
  846.                   SZR = ONE
  847.                   SZI = ZERO
  848.                ELSE
  849.                   TEMPR = C22R / CZ
  850.                   TEMPI = C22I / CZ
  851.                   T = DLAPY2( CZ, C21 )
  852.                   CZ = CZ / T
  853.                   SZR = -C21*TEMPR / T
  854.                   SZI = C21*TEMPI / T
  855.                END IF
  856.             END IF
  857. *
  858. *           Compute Givens rotation on left
  859. *
  860. *           (  CQ   SQ )
  861. *           (  __      )  A or B
  862. *           ( -SQ   CQ )
  863. *
  864.             AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 )
  865.             BN = ABS( B11 ) + ABS( B22 )
  866.             WABS = ABS( WR ) + ABS( WI )
  867.             IF( S1*AN.GT.WABS*BN ) THEN
  868.                CQ = CZ*B11
  869.                SQR = SZR*B22
  870.                SQI = -SZI*B22
  871.             ELSE
  872.                A1R = CZ*A11 + SZR*A12
  873.                A1I = SZI*A12
  874.                A2R = CZ*A21 + SZR*A22
  875.                A2I = SZI*A22
  876.                CQ = DLAPY2( A1R, A1I )
  877.                IF( CQ.LE.SAFMIN ) THEN
  878.                   CQ = ZERO
  879.                   SQR = ONE
  880.                   SQI = ZERO
  881.                ELSE
  882.                   TEMPR = A1R / CQ
  883.                   TEMPI = A1I / CQ
  884.                   SQR = TEMPR*A2R + TEMPI*A2I
  885.                   SQI = TEMPI*A2R - TEMPR*A2I
  886.                END IF
  887.             END IF
  888.             T = DLAPY3( CQ, SQR, SQI )
  889.             CQ = CQ / T
  890.             SQR = SQR / T
  891.             SQI = SQI / T
  892. *
  893. *           Compute diagonal elements of QBZ
  894. *
  895.             TEMPR = SQR*SZR - SQI*SZI
  896.             TEMPI = SQR*SZI + SQI*SZR
  897.             B1R = CQ*CZ*B11 + TEMPR*B22
  898.             B1I = TEMPI*B22
  899.             B1A = DLAPY2( B1R, B1I )
  900.             B2R = CQ*CZ*B22 + TEMPR*B11
  901.             B2I = -TEMPI*B11
  902.             B2A = DLAPY2( B2R, B2I )
  903. *
  904. *           Normalize so beta > 0, and Im( alpha1 ) > 0
  905. *
  906.             BETA( ILAST-1 ) = B1A
  907.             BETA( ILAST ) = B2A
  908.             ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV
  909.             ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV
  910.             ALPHAR( ILAST ) = ( WR*B2A )*S1INV
  911.             ALPHAI( ILAST ) = -( WI*B2A )*S1INV
  912. *
  913. *           Step 3: Go to next block -- exit if finished.
  914. *
  915.             ILAST = IFIRST - 1
  916.             IF( ILAST.LT.ILO )
  917.      $         GO TO 380
  918. *
  919. *           Reset counters
  920. *
  921.             IITER = 0
  922.             ESHIFT = ZERO
  923.             IF( .NOT.ILSCHR ) THEN
  924.                ILASTM = ILAST
  925.                IF( IFRSTM.GT.ILAST )
  926.      $            IFRSTM = ILO
  927.             END IF
  928.             GO TO 350
  929.          ELSE
  930. *
  931. *           Usual case: 3x3 or larger block, using Francis implicit
  932. *                       double-shift
  933. *
  934. *                                    2
  935. *           Eigenvalue equation is  w  - c w + d = 0,
  936. *
  937. *                                         -1 2        -1
  938. *           so compute 1st column of  (A B  )  - c A B   + d
  939. *           using the formula in QZIT (from EISPACK)
  940. *
  941. *           We assume that the block is at least 3x3
  942. *
  943.             AD11 = ( ASCALE*A( ILAST-1, ILAST-1 ) ) /
  944.      $             ( BSCALE*B( ILAST-1, ILAST-1 ) )
  945.             AD21 = ( ASCALE*A( ILAST, ILAST-1 ) ) /
  946.      $             ( BSCALE*B( ILAST-1, ILAST-1 ) )
  947.             AD12 = ( ASCALE*A( ILAST-1, ILAST ) ) /
  948.      $             ( BSCALE*B( ILAST, ILAST ) )
  949.             AD22 = ( ASCALE*A( ILAST, ILAST ) ) /
  950.      $             ( BSCALE*B( ILAST, ILAST ) )
  951.             U12 = B( ILAST-1, ILAST ) / B( ILAST, ILAST )
  952.             AD11L = ( ASCALE*A( IFIRST, IFIRST ) ) /
  953.      $              ( BSCALE*B( IFIRST, IFIRST ) )
  954.             AD21L = ( ASCALE*A( IFIRST+1, IFIRST ) ) /
  955.      $              ( BSCALE*B( IFIRST, IFIRST ) )
  956.             AD12L = ( ASCALE*A( IFIRST, IFIRST+1 ) ) /
  957.      $              ( BSCALE*B( IFIRST+1, IFIRST+1 ) )
  958.             AD22L = ( ASCALE*A( IFIRST+1, IFIRST+1 ) ) /
  959.      $              ( BSCALE*B( IFIRST+1, IFIRST+1 ) )
  960.             AD32L = ( ASCALE*A( IFIRST+2, IFIRST+1 ) ) /
  961.      $              ( BSCALE*B( IFIRST+1, IFIRST+1 ) )
  962.             U12L = B( IFIRST, IFIRST+1 ) / B( IFIRST+1, IFIRST+1 )
  963. *
  964.             V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 +
  965.      $               AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L
  966.             V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )-
  967.      $               ( AD22-AD11L )+AD21*U12 )*AD21L
  968.             V( 3 ) = AD32L*AD21L
  969. *
  970.             ISTART = IFIRST
  971. *
  972.             CALL DLARFG( 3, V( 1 ), V( 2 ), 1, TAU )
  973.             V( 1 ) = ONE
  974. *
  975. *           Sweep
  976. *
  977.             DO 290 J = ISTART, ILAST - 2
  978. *
  979. *              All but last elements: use 3x3 Householder transforms.
  980. *
  981. *              Zero (j-1)st column of A
  982. *
  983.                IF( J.GT.ISTART ) THEN
  984.                   V( 1 ) = A( J, J-1 )
  985.                   V( 2 ) = A( J+1, J-1 )
  986.                   V( 3 ) = A( J+2, J-1 )
  987. *
  988.                   CALL DLARFG( 3, A( J, J-1 ), V( 2 ), 1, TAU )
  989.                   V( 1 ) = ONE
  990.                   A( J+1, J-1 ) = ZERO
  991.                   A( J+2, J-1 ) = ZERO
  992.                END IF
  993. *
  994.                DO 230 JC = J, ILASTM
  995.                   TEMP = TAU*( A( J, JC )+V( 2 )*A( J+1, JC )+V( 3 )*
  996.      $                   A( J+2, JC ) )
  997.                   A( J, JC ) = A( J, JC ) - TEMP
  998.                   A( J+1, JC ) = A( J+1, JC ) - TEMP*V( 2 )
  999.                   A( J+2, JC ) = A( J+2, JC ) - TEMP*V( 3 )
  1000.                   TEMP2 = TAU*( B( J, JC )+V( 2 )*B( J+1, JC )+V( 3 )*
  1001.      $                    B( J+2, JC ) )
  1002.                   B( J, JC ) = B( J, JC ) - TEMP2
  1003.                   B( J+1, JC ) = B( J+1, JC ) - TEMP2*V( 2 )
  1004.                   B( J+2, JC ) = B( J+2, JC ) - TEMP2*V( 3 )
  1005.   230          CONTINUE
  1006.                IF( ILQ ) THEN
  1007.                   DO 240 JR = 1, N
  1008.                      TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )*
  1009.      $                      Q( JR, J+2 ) )
  1010.                      Q( JR, J ) = Q( JR, J ) - TEMP
  1011.                      Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 )
  1012.                      Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 )
  1013.   240             CONTINUE
  1014.                END IF
  1015. *
  1016. *              Zero j-th column of B (see DLAGBC for details)
  1017. *
  1018. *              Swap rows to pivot
  1019. *
  1020.                ILPIVT = .FALSE.
  1021.                TEMP = MAX( ABS( B( J+1, J+1 ) ), ABS( B( J+1, J+2 ) ) )
  1022.                TEMP2 = MAX( ABS( B( J+2, J+1 ) ), ABS( B( J+2, J+2 ) ) )
  1023.                IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN
  1024.                   SCALE = ZERO
  1025.                   U1 = ONE
  1026.                   U2 = ZERO
  1027.                   GO TO 250
  1028.                ELSE IF( TEMP.GE.TEMP2 ) THEN
  1029.                   W11 = B( J+1, J+1 )
  1030.                   W21 = B( J+2, J+1 )
  1031.                   W12 = B( J+1, J+2 )
  1032.                   W22 = B( J+2, J+2 )
  1033.                   U1 = B( J+1, J )
  1034.                   U2 = B( J+2, J )
  1035.                ELSE
  1036.                   W21 = B( J+1, J+1 )
  1037.                   W11 = B( J+2, J+1 )
  1038.                   W22 = B( J+1, J+2 )
  1039.                   W12 = B( J+2, J+2 )
  1040.                   U2 = B( J+1, J )
  1041.                   U1 = B( J+2, J )
  1042.                END IF
  1043. *
  1044. *              Swap columns if nec.
  1045. *
  1046.                IF( ABS( W12 ).GT.ABS( W11 ) ) THEN
  1047.                   ILPIVT = .TRUE.
  1048.                   TEMP = W12
  1049.                   TEMP2 = W22
  1050.                   W12 = W11
  1051.                   W22 = W21
  1052.                   W11 = TEMP
  1053.                   W21 = TEMP2
  1054.                END IF
  1055. *
  1056. *              LU-factor
  1057. *
  1058.                TEMP = W21 / W11
  1059.                U2 = U2 - TEMP*U1
  1060.                W22 = W22 - TEMP*W12
  1061.                W21 = ZERO
  1062. *
  1063. *              Compute SCALE
  1064. *
  1065.                SCALE = ONE
  1066.                IF( ABS( W22 ).LT.SAFMIN ) THEN
  1067.                   SCALE = ZERO
  1068.                   U2 = ONE
  1069.                   U1 = -W12 / W11
  1070.                   GO TO 250
  1071.                END IF
  1072.                IF( ABS( W22 ).LT.ABS( U2 ) )
  1073.      $            SCALE = ABS( W22 / U2 )
  1074.                IF( ABS( W11 ).LT.ABS( U1 ) )
  1075.      $            SCALE = MIN( SCALE, ABS( W11 / U1 ) )
  1076. *
  1077. *              Solve
  1078. *
  1079.                U2 = ( SCALE*U2 ) / W22
  1080.                U1 = ( SCALE*U1-W12*U2 ) / W11
  1081. *
  1082.   250          CONTINUE
  1083.                IF( ILPIVT ) THEN
  1084.                   TEMP = U2
  1085.                   U2 = U1
  1086.                   U1 = TEMP
  1087.                END IF
  1088. *
  1089. *              Compute Householder Vector
  1090. *
  1091.                T = SQRT( SCALE**2+U1**2+U2**2 )
  1092.                TAU = ONE + SCALE / T
  1093.                VS = -ONE / ( SCALE+T )
  1094.                V( 1 ) = ONE
  1095.                V( 2 ) = VS*U1
  1096.                V( 3 ) = VS*U2
  1097. *
  1098. *              Apply transformations from the right.
  1099. *
  1100.                DO 260 JR = IFRSTM, MIN( J+3, ILAST )
  1101.                   TEMP = TAU*( A( JR, J )+V( 2 )*A( JR, J+1 )+V( 3 )*
  1102.      $                   A( JR, J+2 ) )
  1103.                   A( JR, J ) = A( JR, J ) - TEMP
  1104.                   A( JR, J+1 ) = A( JR, J+1 ) - TEMP*V( 2 )
  1105.                   A( JR, J+2 ) = A( JR, J+2 ) - TEMP*V( 3 )
  1106.   260          CONTINUE
  1107.                DO 270 JR = IFRSTM, J + 2
  1108.                   TEMP = TAU*( B( JR, J )+V( 2 )*B( JR, J+1 )+V( 3 )*
  1109.      $                   B( JR, J+2 ) )
  1110.                   B( JR, J ) = B( JR, J ) - TEMP
  1111.                   B( JR, J+1 ) = B( JR, J+1 ) - TEMP*V( 2 )
  1112.                   B( JR, J+2 ) = B( JR, J+2 ) - TEMP*V( 3 )
  1113.   270          CONTINUE
  1114.                IF( ILZ ) THEN
  1115.                   DO 280 JR = 1, N
  1116.                      TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )*
  1117.      $                      Z( JR, J+2 ) )
  1118.                      Z( JR, J ) = Z( JR, J ) - TEMP
  1119.                      Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 )
  1120.                      Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 )
  1121.   280             CONTINUE
  1122.                END IF
  1123.                B( J+1, J ) = ZERO
  1124.                B( J+2, J ) = ZERO
  1125.   290       CONTINUE
  1126. *
  1127. *           Last elements: Use Givens rotations
  1128. *
  1129. *           Rotations from the left
  1130. *
  1131.             J = ILAST - 1
  1132.             TEMP = A( J, J-1 )
  1133.             CALL DLARTG( TEMP, A( J+1, J-1 ), C, S, A( J, J-1 ) )
  1134.             A( J+1, J-1 ) = ZERO
  1135. *
  1136.             DO 300 JC = J, ILASTM
  1137.                TEMP = C*A( J, JC ) + S*A( J+1, JC )
  1138.                A( J+1, JC ) = -S*A( J, JC ) + C*A( J+1, JC )
  1139.                A( J, JC ) = TEMP
  1140.                TEMP2 = C*B( J, JC ) + S*B( J+1, JC )
  1141.                B( J+1, JC ) = -S*B( J, JC ) + C*B( J+1, JC )
  1142.                B( J, JC ) = TEMP2
  1143.   300       CONTINUE
  1144.             IF( ILQ ) THEN
  1145.                DO 310 JR = 1, N
  1146.                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
  1147.                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
  1148.                   Q( JR, J ) = TEMP
  1149.   310          CONTINUE
  1150.             END IF
  1151. *
  1152. *           Rotations from the right.
  1153. *
  1154.             TEMP = B( J+1, J+1 )
  1155.             CALL DLARTG( TEMP, B( J+1, J ), C, S, B( J+1, J+1 ) )
  1156.             B( J+1, J ) = ZERO
  1157. *
  1158.             DO 320 JR = IFRSTM, ILAST
  1159.                TEMP = C*A( JR, J+1 ) + S*A( JR, J )
  1160.                A( JR, J ) = -S*A( JR, J+1 ) + C*A( JR, J )
  1161.                A( JR, J+1 ) = TEMP
  1162.   320       CONTINUE
  1163.             DO 330 JR = IFRSTM, ILAST - 1
  1164.                TEMP = C*B( JR, J+1 ) + S*B( JR, J )
  1165.                B( JR, J ) = -S*B( JR, J+1 ) + C*B( JR, J )
  1166.                B( JR, J+1 ) = TEMP
  1167.   330       CONTINUE
  1168.             IF( ILZ ) THEN
  1169.                DO 340 JR = 1, N
  1170.                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
  1171.                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
  1172.                   Z( JR, J+1 ) = TEMP
  1173.   340          CONTINUE
  1174.             END IF
  1175. *
  1176. *           End of Double-Shift code
  1177. *
  1178.          END IF
  1179. *
  1180.          GO TO 350
  1181. *
  1182. *        End of iteration loop
  1183. *
  1184.   350    CONTINUE
  1185.   360 CONTINUE
  1186. *
  1187. *     Drop-through = non-convergence
  1188. *
  1189.   370 CONTINUE
  1190.       INFO = ILAST
  1191.       GO TO 420
  1192. *
  1193. *     Successful completion of all QZ steps
  1194. *
  1195.   380 CONTINUE
  1196. *
  1197. *     Set Eigenvalues 1:ILO-1
  1198. *
  1199.       DO 410 J = 1, ILO - 1
  1200.          IF( B( J, J ).LT.ZERO ) THEN
  1201.             IF( ILSCHR ) THEN
  1202.                DO 390 JR = 1, J
  1203.                   A( JR, J ) = -A( JR, J )
  1204.                   B( JR, J ) = -B( JR, J )
  1205.   390          CONTINUE
  1206.             ELSE
  1207.                A( J, J ) = -A( J, J )
  1208.                B( J, J ) = -B( J, J )
  1209.             END IF
  1210.             IF( ILZ ) THEN
  1211.                DO 400 JR = 1, N
  1212.                   Z( JR, J ) = -Z( JR, J )
  1213.   400          CONTINUE
  1214.             END IF
  1215.          END IF
  1216.          ALPHAR( J ) = A( J, J )
  1217.          ALPHAI( J ) = ZERO
  1218.          BETA( J ) = B( J, J )
  1219.   410 CONTINUE
  1220. *
  1221. *     Normal Termination
  1222. *
  1223.       INFO = 0
  1224. *
  1225. *     Exit (other than argument error) -- return optimal workspace size
  1226. *
  1227.   420 CONTINUE
  1228.       WORK( 1 ) = DBLE( N )
  1229.       RETURN
  1230. *
  1231. *     End of DHGEQZ
  1232. *
  1233.       END
  1234.